HMOX1 Attenuates the Sensitivity of Hepatocellular Carcinoma Cells to Sorafenib via Modulating the Expression of ABC Transporters

Hepatocellular carcinoma (HCC) represents a common malignancy, and mechanisms of acquired sorafenib resistance during the treatment of HCC patients remain elusive. The present study performed integrated bioinformatics analysis and explored the potential action of heme oxygenase 1 (HMOX1) in sorafenib-resistant HCC cells. Differentially expressed genes (DEGs) of the sorafenib-resistant group as compared to the sorafenib-sensitive group from GSE140202 and GSE143233 were extracted. Fifty common DEGs between GSE140202 and GSE143233 were extracted. Ten hub genes were identified from the protein-protein interaction network based on common DEGs. Experimental results revealed the upregulation of HMOX1 in sorafenib-resistant HCC cells. HMOX1 silence promoted the sensitivity to sorafenib in sorafenib-resistant HCC cells; overexpression of HMOX1 attenuated the sensitivity. In addition, HMOX1 silence downregulated the mRNA expression of ABC transporters in sorafenib-resistant HCC cells, while HMOX1 overexpression upregulated mRNA expression of ABC transporter expression in HCC cells. Further analysis also revealed that high expression of HMOX1 was associated with shorter OS and DSS in HCC patients. In conclusion, our analysis identified ten hub genes associated with sorafenib resistance in HCC. Further validation studies demonstrated that HMOX1 promoted sorafenib resistance of HCC cells via modulating ABC transporter expression.


Introduction
Hepatocellular carcinoma (HCC) represents a common tumor malignancy, and HCC incidence is expected to increase annually [1,2]. Among all the cancer types, HCC is one of the most frequently diagnosed malignancies around the world [1,2]. The HCC patients' overall survival (OS) remains low due to insufficient early diagnosis and the recurrence and metastasis of advanced HCC [3][4][5]. In terms of early diagnosis of HCC, imaging has limitations in diagnostic accuracy and sensitivity, while common serum markers show poor diagnostic performance [6]. With the advent of high-throughput sequencing technology, genetic biomarkers such as cell-free DNA (cfDNA), epigenetic changes, and circulating RNA (microRNAs (miRNAs), long noncoding RNAs (lncRNA), and circular RNAs (circRNA)) from peripheral blood have become the focus of the early diagnosis of HCC [7]. However, there are still limitations in the early diagnosis of HCC. The main limitation of most studies using nucleic acid molecules as biomarkers of HCC is the limited size of the cohort, and the heterogeneity of HCC should be also considered. Sorafenib was first approved by the FDA for HCC therapy at advanced stages. Sorafenib treatment could significantly improve the OS of patients with HCC by 2-3 months. Unfortunately, many patients with HCC had a poor response to sorafenib or exhibited sorafenib resistance after prolonged use of sorafenib [8][9][10]. The molecular mechanisms underlying the acquired sorafenib resistance are complex which may involve epithelial-mesenchymal transition, tumor microenvironment, autophagy, and cancer stem cells [10][11][12].
Besides, studies also proposed that acquired resistance to sorafenib may be associated with the dysregulated signaling pathways including JAK/STAT, PI3K/Akt, and TNFα/NF-κB [10][11][12]. Thus, it is urgent to fully decipher mechanisms associated with sorafenib resistance, which may provide a better therapeutic strategy for treating advanced HCC.  With the rapid development and progress in the highthroughput technologies, the investigation of tumor biology has been focusing on the genomic scale. RNA sequencing has been widely used in identifying novel targets in cancer pathophysiology, including HCC. For instance, Weng et al. performed the RNA-sequencing analysis in HCC and sorafenib-resistant HCC tissues and identified novel cric-FOXM1 as crucial modulator of sorafenib resistance in HCC cells [13]. Wu et al. also carried out global transcriptomic sequencing in sorafenib-resistant HCC cells and emphasized the importance of circRNAs in mediating sorafenib resistance [14]. Wu et al. performed the RNAsequencing analysis and revealed that mitophagy promoted sorafenib resistance through hypoxia-inducible ATAD3A dependent axis [15]. RNA-sequencing studies also proposed that epigenetically activated ADAMTSL5 is a key player in HCC drug resistance [16]. In addition, the integrated bioinformatics analysis of publicly available microarray datasets also provided another strategy for scientists to discover novel mediators associated with the acquired sorafenib resistance in HCC. Liu et al. analyzed dataset GSE109211 and proposed that HCC sorafenib resistance was correlated with identified hub genes and pathways [17]. Jiang et al. analyzed GSE62813, GSE73571, GSE151412, and GSE140202 and found other key mediators in HCC sorafenib resistance [18].
This study employed a strategy to integrate two datasets (GSE140202 and GSE14322) from the GEO database and extracted the DEGs between sorafenib-sensitive and sorafenib-resistant groups. The extracted DEGs were processed for PPI network construction and functional analysis. Furthermore, HMOX1, one of the hub genes, was further validated in sorafenib-resistant HCC cells, and the potential action of HMOX1 in sorafenib-resistant HCC cells was also preliminarily explored.

Materials and Methods
2.1. Extraction of GEO RNA-Sequencing Datasets. RNAsequencing datasets were retrieved from the public GEO database. In the analysis, we searched the RNA-sequencing datasets profiling the mRNA expression between sorafenibsensitive and sorafenib-resistant HCC cells, and samples < 3 in each group were excluded in the analysis. After screening the processed GEO RNA-sequencing datasets using GREIN-iLINCS tool [19], two RNA-sequencing datasets, including GSE140202 and GSE143233, were finally included in our analysis. GPL20795 HiSeq X Ten and GPL16791 Illumina HiSeq 2500 platforms were used in GSE140202 and GSE143233, respectively. In GSE140202, 6 sorafenibsensitive HCC cells (n = 6) and sorafenib-resistant HCC cells (n = 6) were included for analysis. In GSE143233, three samples of HCC tissues and three samples of sorafenib-resistant HCC tissues were included for analysis.

Extraction of DEGs in the Collected RNA-Sequencing
Datasets. The extraction of DEGs in GSE140202 and GSE143233 was performed using the GREIN-iLINCS tool [19]. The DEGs between sorafenib-sensitive and sorafenibresistant HCC cells were extracted in GSE140202; the DEGs between HCC and sorafenib-resistant HCC tissues were extracted for GSE143233. The significant DEGs were extracted using the following criteria: log jFCj > 1:2 and P values < 0.05. The heatmaps of DEGs in GSE140202 and GSE143233 were plotted by using the GREIN-iLINCS tool, and the top 200 DEGs were illustrated in the heatmaps. The R software plotted the volcano plots of the DEGs in GSE140202 and GSE143233 with ggplot function. The Venn diagrams showing the common DEGs between GSE140202 and GSE143233 were also generated using R software.

Functional
Analysis of DEGs between GSE140202 and GSE143233. The GO functional enrichment and KEGG pathway enrichment analyses of common DEGs between GSE140202 and GSE143233 were carried out using EnrichR, an interactive and collaborative HTML5 gene list enrichment analysis tool.

PPI Network Construction Based on Common DEGs
Extracted from GSE140202 and GSE143233. PPI network of the common DEGs between GSE140202 and GSE143233 was generated by STRING [20]. The PPI network was  International Journal of Genomics visualized by using Cytoscape software. To further extract the key submodules in the PPI network, the cytoHubba application from Cytoscape software was applied to visualize submodules derived from the PPI network.       Thermal conditions of amplification on cycles were as follows: 95°C for 2 min, 40 cycles at 95°C for 30 s, 56°C for 30 s, and 72°C for 1 min, then 72°C for 5 min. GAPDH was used as an internal control. Relative expression of the corresponding mRNA was calculated using the 2 -ΔΔCt method. The primer sequences are shown in Supplementary Table 2.

Survival Analysis of Patients with HCC.
The impact of HMOX1 on OS and disease-specific survival (DSS) of HCC patients was evaluated by the Kaplan Meier plotter [21]. Log-rank P < 0:05 was statistically significant.
2.11. Statistical Analysis. The statistical analysis was carried out by GraphPad Prism 5 software (GraphPad Software). All the data were displayed as mean ± standard deviation. Unpaired Student's t-test determined significant differences between treatment groups. P < 0:05 was statistically significant.

Analysis of DEGs between Sorafenib-Sensitive and Sorafenib-Resistant Groups in GSE140202 and GSE143233
Datasets. The DEGs between two groups in GSE140202 and GSE143233 were plotted as heatmaps (Figures 1(a) and 1(b)) and volcano plots (Figures 1(c) and 1(d)). In GSE140202, 542 upregulated and 651 downregulated DEGs were extracted between two groups in GSE140202; in GSE143233, 372 upregulated DEGs and 430 downregulated DEGs were extracted between two groups in GSE143233. The Venn diagram showed that 22 upregulated common DEGs were found between GSE140202 and GSE143233 (Figure 2(a)); 26 downregulated common DEGs were found between GSE140202 and GSE143233 (Figure 2(b)).

Functional Analysis of Common DEGs.
Common DEGs in GSE140202 and GSE143233 were further processed for functional enrichment analysis. In the GO_biological process, common DEGs were enriched in GO terms such as "protection from natural killer cell mediated cytotoxicity"; "antigen processing and presentation of endogenous peptide antigen via MHC class I via ER pathway"; "antigen processing and presentation of endogenous peptide antigen via MHC class I via ER pathway, TAP-independent"; "antigen processing and presentation of exogenous peptide antigen via MHC class I, TAP-independent"; and "myelin assembly" (Figure 3(a)). In the GO_cellular component, common DEGs were enriched in GO terms such as "MHC class I protein complex," "COPII-coated ER to Golgi transport vesicle," "MHC protein complex," and "lumenal side of endoplasmic reticulum membrane" (Figure 3(b)). In GO_ molecular function, common DEGs were significantly  International Journal of Genomics enriched in the GO terms such as "metal ion binding," "serine-type endopeptidase inhibitor activity," "C3HC4type RING finger domain binding," and "apolipoprotein receptor binding" (Figure 3(c)). In KEGG pathway analysis, common DEGs were enriched in "Antigen processing and presentation," "Allograft rejection," "Natural killer cell mediated cytotoxicity," and "Graft-versus-host disease" pathways ( Figure 3(d)).

PPI Network Analysis from Common DEGs in
GSE140202 and GSE143233. PPI network was established based on 48 common DEGs. The 12 nodes and 18 edges were detected in constructed PPI network (Figure 4(a)). Furthermore, the cytoHubba analysis extracted four submodules. The most significant submodule contains four nodes and four edges (Figure 4(b)).
To further explore mechanisms associated with HMOX1-mediated sorafenib resistance of HCC cells, we examined actions of HMOX1 overexpression or silence on ABC transporter mRNA expression in HCC cells. The ABC transporter mRNA expression levels of ABCA6, ABCB1, ABCC1, and ABCG2 in these cells were upregulated 9 International Journal of Genomics compared to their corresponding parental HCC cells (Figure 7). The loss-of-function results demonstrated that HMOX1 silence reduced expression levels of ABCA6, ABCB1, ABCC1, and ABCG2 in these cells (Figure 8). Consistently, the gain-of-function results revealed that HMOX1 overexpression elevated mRNA levels of ABCA6, ABCB1, ABCC1, and ABCG2 in Huh7 and HepG2 cells (Figure 9).

The Prognostic Role of HMOX1 Expression in HCC
Patients. The effects of HMOX1 expression on the prognosis of HCC patients were determined by an online KM-plotter tool, and high expression of HMOX1 was associated with shorter OS and DSS in HCC patients ( Figure 10).

Discussion
The acquired sorafenib resistance has significantly limited the therapeutic potential of sorafenib in advanced HCC, while mechanisms associated with acquired sorafenib resistance remain to be further clarified [10]. In this study, two RNA-sequencing datasets (GSE140202 and GSE143233) from GEO were downloaded for analysis. The DEGs between sorafenib-sensitive and the sorafenib-resistant group were extracted in these two datasets, and a total of 48 common DEGs between GSE140202 and GSE143233 were extracted. Functional enrichment analysis revealed that the common DEGs might be associated with "antigen processing and presentation," "natural killer cell mediated cytotoxicity," "cell adhesion molecules," and so on. Ten hub genes were identified from the protein-protein interaction network based on common DEGs. Experimental results revealed the upregulation of HMOX1 in sorafenib-resistant HCC cells. HMOX1 silence promoted the sensitivity to sorafenib in sorafenib-resistant HCC cells; overexpression of HMOX1 attenuated the sensitivity. In addition, HMOX1 silence downregulated the mRNA expression of ABC transporters in sorafenib-resistant HCC cells, while HMOX1 overexpression upregulated mRNA expression of ABC transporter expression in HCC cells. Further analysis also revealed that high expression of HMOX1 was associated with shorter OS and DSS in HCC patients. Collectively, our results demonstrated that HMXO1 might be associated with sorafenib resistance in HCC.
Analysis of the microarray datasets has become a powerful tool for exploring novel genes that may be associated with cancer progression. In the GSE140202 dataset, TCNOS_00284048 and TCONS_00006019 were highly expressed in the sorafenib-resistant HCC cells compared with parental HCC cells. Knockdown of TCNOS_00284048 and TCONS_00006019 promoted sorafenib-sensitivity of Huh7-SR and HepG2-SR cells [22]. In GSE143233, Lin et al. demonstrated that METTL3 was underexpressed in human sorafenib-resistant HCC and revealed that RNA m 6 A methylation mediated sorafenib resistance via FOXO3-mediated autophagy [23]. Jiang et al. carried out integrated transcriptomic sequencing of GSE140202 and other datasets and identified 13 hub genes and seven promising therapeutic agents for HCC [18]. Li et al. also performed bioinformatics analysis for GSE140202 and found that GINS1 was highly expressed in sorafenib-resistant HCC cells [24]. In our results, we identified 50 common DEGs between GSE140202 and GSE143233, and these DEGs   International Journal of Genomics may be "antigen processing and presentation," "natural killer cell mediated cytotoxicity," "cell adhesion molecules," and so on. Ten hub genes were detected from the PPI network. Among these hub genes, HMOX1 was highly expressed in sorafenib-resistant HCC cells as determined by qRT-PCR assay. HMOX1 is well known for its enzymatic role in regulating cellular homeostasis under stress [25]. Increasing evidence has demonstrated the regulatory role of HMOX1 in cancer progression. For example, Yim et al. showed that HMOX1 is a prognostic marker for bladder cancer recurrence and progression [26]. Park et al. showed that the HMOX1/carbon monoxide axis inhibited transforming growth factor-β1-induced growth inhibition in HCC cells [27]. Inhibiting HMOX1 expression could retard HCC progression via distinct mechanisms [28][29][30]. However, the action of HMXO1 in HCC sorafenib resistance has not been reported yet. Gao et al. showed that inhibition of HMXO1 could sensitize clear-cell renal cell carcinoma to sorafenib [31]. Our study showed that HMOX1 was upregulated in the sorafenib-resistant HCC cells. HMOX1 could promote the resistance of HCC cells to sorafenib. Multidrug resistance is closely regulated by ABC transporters [32]. Studies have demonstrated that sorafenib can interact with ABC transporters such as ABCB1, ABCC1, ABCG2, and ABCC10 [33,34]. In our study, sorafenib-resistant HCC cells exhibited high expression of ABC transporters, including ABCA6, ABCB1, ABCC1, and ABCG2 HCC cells. In addition, HMOX1 silence downregulated the mRNA expression of ABCA6, ABCB1, ABCC1, and ABCG2 in sorafenibresistant HCC cells, while HMOX1 overexpression upregulated these transporters' expression in HCC cells. These results indicated that HMOX1-mediated sorafenib resistance might be associated with the modulation of the expression of ABC transporters.

Conclusions
In conclusion, the present study identified ten hub genes linked to the sorafenib resistance in HCC according to bioinformatics analysis. Further validation studies demonstrated that HMOX1 promoted the sorafenib resistance of HCC cells via modulating ABC transporter expression. However, the bioinformatic analysis was limited to two RNA-sequencing datasets, and future studies should explore more relevant datasets to identify more novel potential mediators in the regulating sorafenib sensitivity of HCC cells.

Data Availability
All the data are available from the corresponding author.

Conflicts of Interest
The authors declare that they have no conflicts of interest.